%% process_fespace
function fespace = process_fespace(mesh, basis_type)
[fespace.node, fespace.elem] = process_refine(mesh, basis_type);
fespace.basis_type = basis_type;
fespace.Nb = size(fespace.node,1);
fespace.Nlb = size(fespace.elem,2);
end

%% process_refine
function [new_node, new_elem] = process_refine(mesh, basis_type)
node = mesh.node;
elem = mesh.elem;
elem_ed = mesh.elem_ed;
edge = mesh.edge;
Np = mesh.Np;
Ne = mesh.Ne;
switch basis_type
    case "P1"
        new_node = node;
        new_elem = elem;
    case "P1b"
        center = zeros(Ne,2);
        for e = 1:Ne
            V = [node(elem(e,:),:)];
            center(e,:) = mean(V);
        end
        new_node = [node; center];
        new_col = (Np+1:Np+Ne)';
        new_elem = [elem, new_col];
    case "P2"
        midpoint = 0.5*(node(edge(:,1),:) + node(edge(:,2),:));
        new_node = [node; midpoint];
        e1 = abs(elem_ed)+Np;
        new_elem = [elem, e1];
    case "P3"
        differences = node(edge(:,2),:) - node(edge(:,1),:);
        point1_3 = node(edge(:,1),:) + differences/3;
        point2_3 = node(edge(:,1),:) + 2*differences/3;
        point_edge = [point1_3; point2_3];
        point_in = zeros(Ne,2);
        mod = [1, 1, 1]/3;
        for e = 1:Ne
            V = [node(elem(e,:),:)];
            point_in(e,:) = mod*V;
        end
        new_node = [node; point_edge; point_in];
        e1 = abs(elem_ed)+Np;
        e2 = abs(elem_ed)+max(e1(:));
        gap_N = max(e2(:));
        idx = elem_ed < 0;
        [e1(idx), e2(idx)] = deal(e2(idx), e1(idx));
        new_col = (gap_N+1:gap_N+Ne)';
        new_elem = [elem, e1, e2, new_col];
    case "P4"
        differences = node(edge(:,2),:) - node(edge(:,1),:);
        point1_4 = node(edge(:,1),:) + differences/4;
        point2_4 = node(edge(:,1),:) + differences/2;
        point3_4 = node(edge(:,1),:) + 3*differences/4;
        point_edge = [point1_4; point2_4; point3_4];
        point_in = zeros(3*Ne,2);
        mod = [2, 1, 1;
            1, 2, 1;
            1, 1, 2]/4;
        for e = 1:Ne
            V = [node(elem(e,:),:)];
            in = mod*V;
            point_in(e,:) = in(1,:);
            point_in(e+Ne,:) = in(2,:);
            point_in(e+2*Ne,:) = in(3,:);
        end
        new_node = [node; point_edge; point_in];
        e1 = abs(elem_ed)+Np;
        e2 = abs(elem_ed)+max(e1(:));
        e3 = abs(elem_ed)+max(e2(:));
        gap_N = max(e3(:));
        idx = elem_ed < 0;
        [e1(idx), e3(idx)] = deal(e3(idx), e1(idx));
        new_cols = reshape(gap_N+1:gap_N+3*Ne, Ne, 3);
        new_elem = [elem, e1, e2, e3, new_cols];
    case "P5"
        differences = node(edge(:,2),:) - node(edge(:,1),:);
        point1_5 = node(edge(:,1),:) + differences/5;
        point2_5 = node(edge(:,1),:) + 2*differences/5;
        point3_5 = node(edge(:,1),:) + 3*differences/5;
        point4_5 = node(edge(:,1),:) + 4*differences/5;
        point_edge = [point1_5; point2_5; point3_5; point4_5];
        point_in = zeros(6*Ne,2);
        mod = [3, 1, 1;
            1, 3, 1;
            1, 1, 3;
            2, 2, 1;
            1, 2, 2;
            2, 1, 2]/5;
        for e = 1:Ne
            V = [node(elem(e,:),:)];
            in = mod*V;
            point_in(e,:) = in(1,:);
            point_in(e+Ne,:) = in(2,:);
            point_in(e+2*Ne,:) = in(3,:);
            point_in(e+3*Ne,:) = in(4,:);
            point_in(e+4*Ne,:) = in(5,:);
            point_in(e+5*Ne,:) = in(6,:);
        end
        new_node = [node; point_edge; point_in];
        e1 = abs(elem_ed)+Np;
        e2 = abs(elem_ed)+max(e1(:));
        e3 = abs(elem_ed)+max(e2(:));
        e4 = abs(elem_ed)+max(e3(:));
        gap_N = max(e4(:));
        idx = elem_ed < 0;
        [e1(idx), e4(idx)] = deal(e4(idx), e1(idx));
        [e2(idx), e3(idx)] = deal(e3(idx), e2(idx));
        new_cols = reshape(gap_N+1:gap_N+6*Ne, Ne, 6);
        new_elem = [elem, e1, e2, e3, e4, new_cols];
    otherwise
        error("Invalid basis type.");
end
end